from numpy import *

def montecarlo_pre(funcion, limites, N):
    # Volumen
    V = prod([b - a for a,b in limites])
    # Puntos random de la funcion y funcion^2 y sus promedios
    suma_f = 0
    suma_f2 = 0
    
    for i in range(N):
        # limites = [(a, b), (c, d), ...]
        puntos = [random.uniform(a, b) for a,b in limites]
        
        f = funcion(*puntos)
        suma_f += f
        suma_f2 += f**2
        
    f_media = (1/N)*suma_f
    f2_media = (1/N)*suma_f2
    # Integral y Error
    integral = V * f_media
    error = V * sqrt((f2_media - f_media**2)/N)
    
    return integral, error

def montecarlo(funcion, limites, N):
    
    V = prod([b - a for a,b in limites])
    dimension = len(limites)
    puntos = random.rand(N, dimension)
    
    for i,(a,b) in enumerate(limites):
        puntos[:,i] = a + (b-a)*puntos[:,i]

    f = funcion(*puntos.T)

    f_media = (1/N)*sum(f)
    f2_media = (1/N)*sum(f**2)
    
    integral = V * f_media
    error = V * sqrt((f2_media - f_media**2)/N)
    
    return integral, error

def f(x):
    return  (1 - x**2)**(1.5)

lim = [(0, 1)]

print("N", "I", "Error", sep="                      ")
print("-"*55)
for n in range(2, 9):
    N = 10**n
    I, error = montecarlo(f, lim, N)
    print(f"{N:>12}\t{I:20.12f}\t{error:20.2e}")
print("-"*55)